/*

calculation of covariance and correlation for pairs of vectors in N-D matrices

FORMAT:       [cv, r] = cov_nd(X,Y)

Input fields:

      X, Y        N-D double matrices of same size

Output fields:

      cv          covariance matrix
      r           correlation coefficients matrix

Note: cv and r are computed over the LAST dimension, so if X and Y
are 10-by-30-by-50 matrices, cv and r are 10-by-30 matrices.

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

*/

#include "mex.h"
#include "matrix.h"
#include "math.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	int docr, c, ic, nd, nl, dr[64], rd, rs, ts;
	const int *d1, *d2;
	double dnl, dnx, *p1, *p2, *t1, *t2;
	long double m1, m2, x1, x2, cv;
	double mp;

	if (nrhs != 2 || nlhs > 2)
		mexErrMsgTxt("Bad number of input/output arguments.");

	if (nlhs > 1)
		docr = 1;
	    else	
		docr = 0;

	nd = mxGetNumberOfDimensions(prhs[0]);
        if (nd != mxGetNumberOfDimensions(prhs[1]) || nd > 64)
		mexErrMsgTxt("Matrices must have the same size.");

	d1 = mxGetDimensions(prhs[0]);
	d2 = mxGetDimensions(prhs[1]);
	for(c = 0; c < nd; ++c)
		if (d1[c] != d2[c] || d1[c] == 0)
			mexErrMsgTxt("Matrices must have the same size and not be empty.");

	for(c = 0; c < 2; ++c)
		if (!mxIsNumeric(prhs[c]) ||  mxIsComplex(prhs[c]) ||
		     mxIsSparse(prhs[c])  || !mxIsDouble(prhs[c]))
			mexErrMsgTxt("Matrices must be numeric, real, full and double.");

	rd = nd - 1;
	if (d1[rd] == 1)
		mexErrMsgTxt("The last matrix dimensions must not be singleton.");

	p1 = (double *) mxGetPr(prhs[0]);
	p2 = (double *) mxGetPr(prhs[1]);

	rs = 1;
	for(c = 0; c < rd; ++c) {
		dr[c] = d1[c];
		rs *= dr[c];
	}
	nl = d1[rd];
	ts = rs * nl;
	dnl = (double) nl;
	dnx = dnl - 1.0;

	plhs[0] = mxCreateNumericArray(rd, dr, mxDOUBLE_CLASS, mxREAL);
	t1 = (double *) mxGetPr(plhs[0]);

	if (!docr) {

		for(c = 0; c < rs; ++c) {
			m1 = 0.0;
			m2 = 0.0;
			mp = 0.0;
			for(ic = c; ic < ts; ic += rs) {
				m1 += p1[ic];
				m2 += p2[ic];
				mp += p1[ic] * p2[ic];
			}
			m1 /= dnl;
			m2 /= dnl;
			mp /= dnl;
			*t1++ = (double) ((mp - m1 * m2) * (dnl / dnx));
		}

	} else {

		plhs[1] = mxCreateNumericArray(rd, dr, mxDOUBLE_CLASS, mxREAL);
		t2 = (double *) mxGetPr(plhs[1]);
		for(c = 0;c < rs; ++c) {
			m1 = 0.0;
			m2 = 0.0;
			x1 = 0.0;
			x2 = 0.0;
			mp = 0.0;
			for(ic = c; ic < ts;ic += rs) {
				m1 += p1[ic];
				m2 += p2[ic];
				x1 += p1[ic] * p1[ic];
				x2 += p2[ic] * p2[ic];
				mp += p1[ic] * p2[ic];
			}
			m1 /= dnl;
			m2 /= dnl;
			x1 /= dnl;
			x2 /= dnl;
			mp /= dnl;
			cv = (mp - m1 * m2) * (dnl / dnx);
			*t1++ = (double) cv;
            cv /= sqrt((x1 - m1 * m1) * (x2 - m2 * m2)) * (dnl / dnx);
            if (cv < -1.0)
                cv = -1.0;
            else if (cv > 1.0)
                cv = 1.0;
			*t2++ = (double) cv;
		}
	}
}
